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A theoretical analysis [Angelani et al., Phys. Rev. Lett. 96, 065702 (2006)] predicts glassy be- 
CIh' havior of light in a nonlinear random medium. This implies slow dynamics related to the presence 

of many metastable states. We consider very general equations (that also apply to other systems, 
like Bose-Condensed gases) describing light in a disordered non-linear medium and through some 
^sj ' approximations we relate them to a mean-field spin-glass-like model. The model is solved by the 

replica method, and replica-symmetry breaking phase transition is predicted. The transition de- 
scribes a mode-locking process in which the phases of the modes are locked to random (history 
and sample- dependent) values. An extended discussion of possible experimental implications of our 
analysis is reported. 

I. INTRODUCTION 

In a nutshell, laser action in a stochastic resonator (SR) defines a random laser (RL). Following the original 
Lethokov's article, |l| a SR is a disordered medium sustaining a large number of electromagnetic modes with over- 
lapping resonances. The modes are not necessarily localized (in the Anderson sense), but can be extended modes in 
a random medium; they have typically a finite life-time and are sometimes referred to as "quasi-modes" . Generally 
speaking, we will refer to a RL as a multi-mode laser system that displays some disorder; this will be described by 
a probability distribution and we will then consider different realizations of the system. Such a general definition 
not only embraces the experiments addressed below, but also include standard lasers with a disordered cavity, or 
integrated devices, as for example ordered photonic crystals Q] infiltrated by some active (i.e. doped) soft-material, 
like liquid crystals or polymers, that induces a given amount of disorder, or even intentionally disordered photonic 
J^t . crystals enriched by quantum wells providing optical gain. 

In the early developments, the theoretical framework at the basis of RL has relied on light diffusion pUj-Lj These 
studies stimulated many investigations concer ning photon dynamics in a disordered medium, up to considering the 
I" ■ quantum transport of photons |6ll7ll8ll9l llCl lll. 12]. Subsequent detailed numerical studies revealed how important for 
RLs is the nature and the distribution of localized modes in random amplifying media, in particular in the strongly 
Zr ■ scattering regime jlj, uM- Experiments were reported on the emerging of many coupled oscillation modes while 
increasing the pump energy and the consequent non-trivial dynamics of the resulting optical signals |15l lid , Il7l Il8| . 
Coupling of modes were addressed in [l8J ] , as the fact that the maximum observed number of modes increases with 
the pumping intensity and with the sample volume |l9l |20H . Recent results pointed out new key issues concerning the 
physics of random lasers, as the role of extended modes [2l],|22], or the presence of specific fluctuations. [23. |24| . 

When considered from a semi-classical perspective, a multi-mode random laser strikingly displays those ingredients 
which are typical of the physics of complexity: i.e. randomness and nonlinearity. The latter is due to typical mode- 
interaction processes, like mode-competition and mode-locking J25L 129.1271 128 . |29| . Complex processes in laser physics, 
including nonlinear optics, are well known and studied (see e.g. |3(J,|3l|), up to recent investigations in multi-mode 
systems |32> [33|, [34( and successful reformulations of standard laser thermodynamics 35, 3al37ll38l |39J | . The extension 
of these approaches to RL immediately leads to the application of the statistical theory of disordered systems, of which 
spin glass theory is a paradigm [401 ] . and which is the subject of the present manuscript. 

We show that, in the presence of a large number of coupled modes (extended or not), a mode-locking (ML) process 
can be observed. ML is related to the relative phases between resonant states, which in some cases become locked at 
the same value. For a standard laser it can be realized by an active device, like an acusto-optic modulator, or can be 
self-starting as in the presence of a nonlinearly mediated mode interaction J4lj] . We can expect that RL-ML appears 
when many modes are put into oscillations, and their amplitudes are clamped at the oscillation values, which are 
random variables. The temporal dynamics of the emitted signal is indeed strongly related to the phases [421 ] . given 
the fact that the mode amplitudes vary on a much longer time scale. The latter circumstance favors the consideration 
of the mode amplitudes as "quenched" (i.e. random but slowly varying) variables, and the phases are to be taken as 
the relevant dynamical variables. The mode-locking process in standard lasers is now recognized as a thermodynamic 
phase transition [35ll3a . l37ll3a . l39l ]: it is expected, therefore, that the mode- locking transition for a RL takes the form 
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of a phase transition in a disordered system. 

This manuscript follows a recent letter |43|. and furnishes: extensive and new details on the derivation of the 
analytical results (including the stability analysis that was previously not reported); a discussion of the underlying 
working hypotheses; a discussion on the nature of the considered electromagnetic modes; the analysis of possible 
experimental frameworks where glassy behavior of light can be observed. The paper is structured as follows. In 
section II, we will review coupled mode theory in a dielectric resonator in the presence of a nonlinear susceptibility; in 
section III we will specialize the approach by deriving the leading model for our analysis; in section IV we will apply 
the methods used in spin glass theory to solve the model; section V is focused on a discussion of the physical meaning 
of our results, using real-world units, and of possible experimental setups; conclusions are drawn in section VI. 

II. COUPLED MODE THEORY EQUATIONS 

The physical system under consideration is an open electromagnetic cavity supporting modes at optical frequencies. 
The cavity is characterized by the presence of disorder; for example randomly structured dielectrics in between a 
couple of mirrors, or a mirrorless system (e.g. a distribution of dielectric particles) such that there is a sufficiently 
high refractive index contrast, so that localized modes (which means modes belonging to the discrete spectrum of 
the eigenvalue problem given by the Maxwell equations, as detailed below) do exist. This is the case, for example, 
of a disordered distribution of TiC>2 particles, of semiconductor powders in a liquid or a glassy matrix, or of a 
nanostructured microcavity filled by a randomly fluctuating material like liquid crystals or soft matter. The localized 
modes supported by these systems can be very different, depending on the degree of localizations, e.g. they can be 
distributed over all the the dielectric sample (as those investigated for example in the experiments reported in |2l|) or 
correspond to well localized states (as those numerically analysed in [lj]); this distinction can influence the properties 
of the interaction between modes, leading to interesting effects, as will be discussed in the following. However the 
physical picture we will obtain is expected to be independent of the details of the interaction. 

Models for multimode nonlinear optical cavities have been largely reported in literature (see e.g. @>II1jE3)- Typ- 
ically they result into coupled equations for complex amplitudes, which can be obtained using various and equivalent 
approaches. In order to fix the notation and for the benefit of the non-expert reader, here we will briefly report a 
derivation based on a multiple scale approach. The electromagnetic cavity [a dielectric resonator (DR)] , is described 
by a (static) refractive index profile n(r). Such a kind of system may support the existence of resonance modes, which 
can be either localized or distributed in the system. Maxwell's equations are written as 

VxH = e n 2 (r)<9 t E m 

V x E = -Mo<9 t H W 

The electric and magnetic fields can be expanded in normal modes with angular frequencies w n and eigenvalues E„ (r) 
and H n (r) as: 

E = RcE„ E„(r) exp(-iw n f)] ,„. 

H = ReE„ H„(r) exp (-«,;»*)] [Z > 

The latter quantities satisfy the generalized eigenvalue problem 

CT S = losMFs (3) 

while being 
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Given a volume V much wider than the DR, over which periodical (Born- Von Karman) boundary conditions are 
posed, and introducing the complex valued scalar product 

(A,B) = f A* BdV (7) 

Jv 

it turns out that C and M. are self-adjoint operators. As a result uj n are real valued and the eigenvectors are orthogonal 
with weight M.. Furthermore, since (E„, H„) and (E* , — H* ) correspond to the same eigenvalue oj n , E n can be taken 
as real valued. 

The average electromagnetic energy for each un-normalized mode is given by 

E. = ~ / e Q n 2 (r)\E s \ 2 + vo\Ks\ 2 dV = ] (T S ,MT S ) . (8) 



In the following the modes are normalized in such a way 

^(T s ,MT q ) = S sq (9) 

Next we consider the perturbed Maxwell equations in the presence of a nonlinear polarization Pjvx, such that the 
overall dielectric displacement vector is given by D = e n 2 (r)E + Pnl and J = dtPjyL a generalized current: 

V xH = e a n 2 (r)d t E + i 1 J 

V x E = -McAH, [W) 

r\ is a bookkeeping perturbation parameter to be set equal to one at the end of the derivation. Our aim is to write 
the solution of the nonlinear Maxwell equations as a superposition of modes such that the leading order has the form 



E = ReE n v^o n (t)E n (r) exp(-iu n t)] 



(11) 



H = Re E„ v^ a n(*) H n(r) exp(-iuj n t)} 
and the complex amplitudes a s are such that the total energy stored in the DR is 

C 2-i m G m Zj nl 0J m \& m | . V J 

There are various techniques to derive the leading equations for the a s , here we adopt the multiple scale method 
(see e.g. 46]). The perturbative expansion is written as (with obvious notation) 



E = Re{X; r JV^ra n (?7i,77 2 t, ...)E„ + rjE£ ) + ...} cxp(-iw„i)} 
H = Rc{J2 n [V^an{vtiV 2 t, ...)H„ + rjH^ + ...]exp(-iu n t)} 



2+ m , „tj(i) . i ... ( 13 ) 



where the amplitudes are taken to be dependent on the slow scales t n — rj n t, as the first and higher order corrections 

i the temporal derivatives are written as dt — dt + f]dt^ + ....Letting 

Pnl =Re[^P„(ii,t 2 ,-)exp(-MJ n t )] (14) 



like E n , the fastest scale is to = t and the temporal derivatives are written as dt — dt + rjdt 1 + ....Letting 



and 

J = Re[^J„(ii,t 2 ,...)exp(-w n t )] =d t PwL (15) 



with 






(16) 



it is 

j£°) = -iUnPW. (17) 

Using the previous machinery into the nonlinear Maxwell equations at the first order in r\ it is found for the term 
oscillating with exp(— iw s t ) 

OF<P - UJsMfP - B s (18) 



while having 
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and 
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The Fredholm theorem applied to l|18|) . states that the solvability condition is the orthogonality with the kernel 
solution, i.e. J- s : (J- S ,B S ) — 0. Hence 

^ = ^(E S ,P(°)). (21) 



dh 4 

Going back to the original variables, we have the desired result 

da s (t) y/UJ, 



dt U Jv K(r)-P s (r)dV. (22) 

III. NONLINEAR SUSCEPTIBILITY AND MODE INTERACTIONS IN ACTIVE RANDOM CAVITIES 

We consider the case in which many modes are put into oscillations and interact due to the nonlinearity of the 
amplifying medium. In resonant systems the nonlinear optical response can be found from the density matrix equations 
in a two-level system, as originally investigated by Lamb |26( . The component of the nonlinear susceptibility oscillating 
at uo s is modelled as usual [4jj and is written as: 

p s = Yl Xa0 1 s(uJs]u q ,uj r ,-Lj p ,r)EP(r)E^(r)Ef.(r) y /IJ^J^j;a (! a r a* p (23) 

UJ s -\-LJp=UJ q -\-UJ r 

where x is the third order response susceptibility tensor, which in general depends on the positions in the DR. Using 
(I23JI the coupled mode theory equations ill'L'l) read as 

8 = ~7;^29spqra q a r ap , (24) 



dt 2 

pqr 



while being 



g spqr = Vggp^ I XaPlS (cj s ;u; q: u; r ,-cj p ,r)Ef(r)E^r)E^(r)E 5 r (r)dV . (25) 

« Jv 

A. Mode interactions and the role of localized modes 

Our treatment follows early works on multimode cavities [261 |28| and consistently, in the previous equations, 
intermode frequencies and higher harmonics are neglected because they have in general a lower Q-factor if compared 
to those of the supported cavity modes. Additionally, since the sum in Eq. I|24|l is extended to all the modes 
combination satisfying the condition lu s — cu q + tu r — lj p , we recall that the frequencies satisfying this relation can be 
divided into three categories |23|: (a) u> s = ui q and ui r — u pi (b) oj s = u> r and u> q = uu p \ and (c) uu s = u> q + u r — uo p 
excluding (a) and (b). Categories (a) and (b) were shown to determine the oscillation values of the energies of the 
modes £ s , and provide terms like self and cross-saturation, as also recently considered in [2(|, with reference to RLs. 
The third group are the "combination tone terms" |28| which were originally neglected, even if it was later recognized, 
through numerical calculations, to have a role when the number of modes increases |43|. We are interested to the 
regime in which a large number of modes is put into oscillation in a limited spectral range around a given carrier 
wavelength ujq (which can be taken as the resonant angular frequency of the active medium) and we will show below 



that in RLs, the combination tone terms provide a complex structure to the laser dynamics, as due to the fact that 
for an increasing number of modes they couple almost all the cavity resonances. 

In general, the resonant condition for the mode- locking processes ui s = u) q + oj r — oj p , does not need to be satisfied 
exactly but in such a way that the mode combination tone u> q + u) r — lu p lies within the lincwidth at uj s (this is 
discussed e.g. in |29| with reference to three modes mode- locking). In the presence of many modes oscillating in 
a small bandwidth, and such that the linewidths are overlapping, as it is typically the case for RLs (see the cited 
references) and (by definition) for SRs, many mode combination tones will couple to u> s , for which we have taken 
lu s = uj q + uj r — ojp in (|24fl . This opens the way to a "mean field theory" where all the modes are coupled, i.e. the sum 
in (|24() is over all the possible values of pqr. We will describe this regime, moreover considering the thermodynamic 
limit as the number of modes goes to infinity. 

However, it is worth to observe that the coupling g spqr in 1251) is related to the spatial overlap of the four modes 
E s , E p , E q , E r that enter in the integral. In the case of extended modes, all the modes will have large spatial overlap 
and g spq r 7^ for all spqr, so the "mean field limit" above is expected to be a very good approximation. On the 
contrary, in the case of strongly localized modes with localization len gth £, i t is reasonable to expect that only a finite 
number of modes will be supported in a localization volume ~ £ 3 [lUll4l . ll8J |. so that the coupling g sp qr will be nonzero 
only for those modes which are large in the same (or in adjacent) localization volumes. In this case the interaction will 
be short range, i.e. in the sum (|24|l only quadruplets of "nearby" modes will appear. Many intermediate situations 
between the extended and the strongly localized ones might happen in random lasers [II], [lj, Il8l l2ll | and indeed the 
precise nature of the modes in these systems is not completely clear. 

In the short range case, the basic phenomenology of the glass transition we will find (slow dynamics, random mode- 
locking) remains the same, but the physics of the system is strongly affected by activated processes (nucleation, barrier 
crossing, etc.) which are negligible in the mean field limit. Indeed the nature of the glass phase of short range spin 
glasses is still a debated problem [43. Note that the localization length (and thus the interaction range) may vary 
on many orders of magnitude and can be experimentally controlled 11], at variance to what happens in spin glasses 
and molecular or colloidal glasses, where the interaction range is fixed by the property of the material and is always 
of the order of the interparticle distance. This observation opens the way toward the possibility of an experimental 
investigation of the crossover between the mean field limit and the short range case that might be crucial for the 
theoretical understanding of the glass phase in short range systems. 

To summarize, we will assume that i) all the lasing modes have frequency lo s = ujq, luq being the resonant frequency 
of the active medium, so that the constraint u> s — Lu q + uj r — uj p can be released, and that ii) the spatial overlap of the 
modes is large, so that the integral in (|25|1 will be not negligible for any quadruple of modes. Under these hypotheses 
a "mean field" treatment in which all quadruples of modes interact will be a very good approximation. Nevertheless, 
we expect the physical picture we will find in the following to hold under much more general assumptions on the 
interaction between modes. Its modifications due to the violation of the hypotheses above will be very interesting for 
the theory of spin glass systems. 

B. The "quenched" approximation: a Langevin equation for the phases 

Letting a s (t) — A s (t)exp[iip s (t)], we take A s as slowly varying with respect to (f 8 . Indeed, the facts that the 
temporal variation of the phases is on a time scale faster than that of the amplitudes, and that fluctuations in a 
cavity take place because of the random interference between modes and not because of the intensity fluctuations of 
individual modes, are well established from the theor y o f mode-locking of standard multi-mode lasers 1271 142. l44|. 
Previous analytical, or semi-analytical, studies [2g,|23,|28| (the RL case has been recently considered in J2(J]) relayed 
on the so-called "free run approximation" , i.e. the phases are taken to be rapidly varying and independent and can 
be averaged out (see Appendix 0. This turns out into removing all phase-dependent terms in (|24|) and the resulting 
equations determine the amplitudes A Sl and hence the energy into each mode £ s , which stays clamped at this value 
after that the corresponding mode has been put into oscillation. As far as the phases can be taken as independent, 
the output laser signal displays small oscillations around an equilibrium value, because the noises into each mode 
amplitudes are independent. It is clear that in this approximation the combination tone terms in (|24|l will disappear 
due to the averaging over the phases. However, since the beginning J2y, [2£| (and later also confirmed by detailed 
numerical investigations |48j) it has been known that this regime holds as far as beating between modes, due to the 
mode combination tones, are negligible; and this is valid if a few modes with nonoverlapping resonances are excited. 
Conversely, the mode combination terms are known to be responsible of "mode locking" processes, that in standard 
laser provide a fruitful approach to the generation of ultra-short pulses (42, |4J] . 

Gain (described by an amplification coefficient 7 S ) and radiation losses (measured by a s ) are included in the equation 



of motion for the complex amplitudes following a standard approach [44( : 

—TT = -^y^gspgragQrQp + (is ~ a s )a s + T] s (t) 



(26) 



pqr 



having introduced, as usual, a complex noise term, mainly due to spontaneous emission (see e.g. [50t |5T|), with 
(VpWVqit')) = (Vp(i)r]*(t')) = and (r]p(t)r]*(t')) = 2k B T bath 5 pq S{t - if), with k B the Boltzmann constant and T bath 
an effective temperature, whose expression will be reported in a later section. In (|26[) . the sum has been extended 
over all the modes, as discussed above, and the contribution of each possible combination tone to the amplitude a s is 
given by the relevant coupling coefficient g spqr - 

The tensor g is a quantity symmetric with respect to the exchange of s «-> p, q <-> r, while under {s,p} *-* {g, r} 
one has g spq r = g'q rsp , see Eq. (l2*3l and J2o,[^Z|. Introducing the real- valued potential function 



H = -Re 

4 



g r spqr&qQr® p Q' s 



, / j 9spqr a q a r a p a s ,■ / y 9spqr a q a r a p a 



and letting H = X) s ( a s ~ 7s)|a s | 2 + H , the resulting model IJ2l)|) is re- written as 

da s dH. 



where 



dt 



d _ 1 
d^* ~ 2 



(9n 
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da R da 1 



(27) 



(28) 



(29) 



The previous equation can be cast in the form of a standard Langevin equation for a system of N "particles" moving 
in 2 AT dimensions (represented by {a^,a I s \ s =i..N) lilii Lil| and its invariant measure is given by exp(— 'H/ksTbath)- 

The simplest case is attained when g can be taken as real valued. Indeed, considering Lamb theory for a two 
level system J26J, which is the only approach providing an explicit expression for the susceptibility tensor \, on e can 
show that the imaginary part of g vanishes as all the resonant frequencies are packed around a given value loq. The 
generalization to a complex g is discussed in Appendix fXl 

Finally, the phases ip s can be taken as the relevant dynamic variables, due to the quenched approximation for the 
amplitudes A s , see Appendix 1X1 and the Hamiltonian is written as 



H(G, if) = n o + ^2 G spqr cos(if s + ip p -ipq-ip r ) 



(30) 



{sp},{qr} 



where U = J2 S ( C 

*-* srtai £<g unarms <t*-r>s±a<f*- 



7 s )v4g is an irrelevant constant term (as long as the amplitudes A s are constant) and 
spqr—"}j8pqr-n-s<™-p-™-q<™-r is the real-valued coupling. Note that the couplings G spqr are symmetric under internal 
permutations of the sets {s,p} and {g, r} and also under exchange {s,p} <^ {q, r}. Indeed the couplings have the 



{sp},{qr} 



same symmetry of the interaction term cos^s + ip p — ip q — ip r ). To count each term only once, the sum ^2 

(|30|l has been restricted only to the values of spqr which are not related by the symmetries above, and correspondingly 

a factor of 8 has been added in the coupling. 

Hereafter we will consider these G coefficients as "quenched" (due to the slow t dependence of A s ), and the relevant 
phase space is reduced to that spanned by ip s . The pump energy which controls the average energy into each mode 
(and hence the amplitudes A s ) fixes the amplitude of G. 



C. Gaussian random couplings 



If the cavity is realized by a random medium, as described above, the coupling coefficients g are random variables, 
i.e. they will depend on the specific sample one is considering. For a given cavity realization, the values of the coupling 
coefficients g are determined by the specific nonlinear mechanism, i.e. by the function \ m l|23|) . the mode frequencies 
oj s and amplitudes A s and by the mode profiles E s (r) as expressed in (|25f) . All these are sources of randomness in 
the computation of g. However, as we cannot compute the properties of the model for a specific choice of the g, we 
will assume that the couplings are drawn from a given probability distribution and compute average properties of 
the system with respect to this distribution. It is possible to show that, in the thermodynamic limit, these average 



quantities will be representative of many of the properties of a given sample (e.g. the free energy), see the discussion 
in next sections and |4Cj . 

Given the fact that the fields are real-valued functions with positive and negative values at each point in the medium, 
and additionally, parity of the modes may eventually make some coupling vanishing, a possible choice, in order to 
simplify the problem as much as possible, is to consider the signs of these coupling coefficient (those corresponding 
to the mode-combination tones) as random and treat them as Gaussian independent variables with zero mean, as 
detailed below. This choice is further supported by the fact that, in general, the mode frequencies are symmetrically 
distributed with respect to the resonant frequency luq, and correspondingly the sign of the nonlinear susceptibility 
larg ely varies |4jj . The hypothesis of zero mean can be removed by generalizing the treatment reported below following 
|40| , leading to a very rich phase diagram 52j ; different distributions of the couplings can also be investigated but the 
problem becomes more difficult. 

Additionally, it is important to point out the scaling properties of the Hamiltonian 13Ufl . Recalling that 
Ef = 0(V^ 1/>2 ) (due to the normalization) and Xq/37<5 — 0(1) are random variables, one has, from Eq. 125J I. 
g ~ V~ 2 f v R(r)dV ~ V~ 3 ^ 2 , as R(r) = V 2 \E E E E is an 0{\) random variable whose integral scales as V 1 ! 2 . 
The coupling Gs then scale as (A 2 ) 2 goV~ 3 / 2 . By a simple rescaling, the invariant measure can be written as 
exp[— [3H( J, <p)], where J pqrs — G spqr /(go(A 2 ) 2 ) has standard deviation \jV 3 l 2 ex 1/N 3 / 2 , as the number of modes is 
proportional to the volume of the cavity, see e.g. 18J : then, conventionally we will set (J 2 ) = 8/iV 3 and include all 
the system-dependent constants in the definition of /3. Note that this scaling of the Js guarantees that the Hamilto- 
nian is extensive, i.e. the average energy is proportional to volume 40, 53] . The parameter that controls the phase 
transition is then (3 = 1/T = V 2 /kBTt a th, where V 2 = (A 2 ) 2 g . We recall that uj (A 2 ) measures the average energy 
per mode, while go is a material-dependent constant. Then V is proportional to the energy stored on average into 
each mode. Hence, the relevant parameter for the lasers model is the adimensional "temperature" T: lowering T can 
be obtained both lowering the bath temperature Tf, at h (e.g. acting on the noise, as done in recent experiments on the 
thermodynamics of standard lasers [3jj) or increasing "the pumping rate" V . 

IV. REPLICA ANALYSIS OF THE MODEL 

Our interest here is to draw a mean field statistical description for random lasers, which enables to go beyond the 
"free run approximation" and unveil the complex structures of the states of these systems, due to "random mode- 
locking" processes. We can do this by computing the thermodynamic properties of the Hamiltonian (|3U[I describing 
the stationary states of the system. 

Summing up, the random laser studied in the above sections is described by the disordered mean-field Hamiltonian 

H = ^2 Jspqr COs(lfi s +tp p — <p q — <fir) (31) 

{pr},{sp} 

where {y s } are angular variables, ip s € [0, 2tt), and J sp q r are independent Gaussian random variables with zero mean 
and variance J 2 = a 2 = 8/N 3 . The sum in i|31|) is restricted to the values of spqr that are not related by the 
symmetries of the interaction term, see the discussion after IJ30II . Our purpose is to study the thermodynamics of this 
model. In particular we are interested in average properties of the model considering the variables Js as quenched: 
that means, we will average the free energy, and not the partition function, over the distribution of the couplings: 
averaging the partition function correspond to considering the Js as dynamical variables evolving on the same time 
scale of the phases. The average of the free energy can be done by means of the replica trick [I(J]. Here we report the 
calculation in full detail, even if it is very similar to the replica calculation for the p-spin model, see e.g. J53|. 

A. Replicated partition function 

The partition function is 



Z N ((3,J)= d{cp}e-^ J ^, (32) 



with free energy 



f N (J3,J)=~biZ N (/3,J). (33) 



We want to calculate the free energy averaged over the disorder 



[z N {p,j)Y-\ 



-pf(0) = lim -lnZjvOS, J) = lim lim L u ' J , (34) 

N-too TV Af-toon^O 7liV 

where the overbar denotes the average over the random coupling Js and, as usual in the replica method, one uses the 
formula In Z — \im n ^o(Z n — l)/«, introducing the partition function of n independent replicas Z n (J) = [Zn(/3, J)] n 
with the same random couplings Js. It is possible to show J4(j that the free energy is self-averaging, i.e. one has 

lim f N (0,J)=f(0) (35) 

with probability one with respect to the distribution of the couplings Js; in other words, in the thermodynamic limit 
the free energy of a given sample is given, with probability 1, by the average of the free energy over the disorder, that 
we are able to compute. 

In the following we will neglect all the multiplicative constants growing as powers of N in the partition function 
as they do not contribute to the free energy. The Js have distribution P(J) = y / N 3 /l6ncxp (— J 2 iV 3 /16): by the 
relation 



j dJP(J)e AJ = exp 



4** 

N 3 



(36) 



for integer n, one has 



Z n {J) = / [J d ^"> M- ff -"fc?> , (37) 



with 



H *ff=NsYl E cos(^ + ^-^-^)cos(^ + ^-^-^) 

a.b {sp},{qr} 

8/V- 3 1,JV 

" Yl E cos (^ + ( p a p-<- tf) cos (^ +4-$- & (38) 



V2 "/ a.b spqr 

N 



■E [l^l 4 + \ R ab 



[\^ab\ - |-"-a&| J 
b 



where in the second line the constraints on the sets {s,p}, {q, r} have been released. In the above expression we have 
introduced the quantities: 

Qab = N- 1 Y,e i(tp *-' p * ) , (39) 

i 

R ab = N~ l J^ e i( *'" + *'' ) . (40) 

i 

Note that Q aa = 1 and <3& a = Q* a bi while i?f, a = i? a fc. Also note that the effective Hamiltonian H e ff depends only on 
the global variables Q a b and R a b- 

Using the notation S(z) = 5(z R )S(z I ), the partition function can be written as 

z%j) = / n dqab n dr - b e*-*-"^'^ f n ^ a ) n *(«»& - &*) n ^ r ^ - ^ > ^ 

■^ a>fc a>b J a=l a>b a>b 

and the second integral can be easily evaluated introducing the integral representation for the complex <5-function 

S( z ) = f JL^ e ^{-n , dl = d^dl 1 (42) 

J (27r)^ 

and the integral is done on the imaginary axis of the complex l R , I 1 planes (i.e. one has to consider both l R and I 1 
as complex numbers). With some algebra one gets, with the convention that [ab] — > a > b and (ab) — > a > b, and 



summing over the repeated indexes, 

/ d{<P a } Yl S ( q ( ab ) ~Q(ab))Y[ S ( r [ab] ~ R[ab]) = 



where 



(ab) [ab] 

dl(ab)dn[ab] exp iVRe(Z* ab) g (a6) + n\ab\ r [ab]) + N\nZ(l( ab) , /.i [ab] ] 



-Tie (if , P i ( i P a ~ i P b ) _|_ n* ,p i (v a +V b ) 
lxe y(ab) e +M[a6] e 



(43) 



Z{l(ab)^[ab\) = I d[ip a ]exp 



(44) 



(note the difference between d{(p a } — J\ a * dipf and d[ip a ] — J\ a d(p a ). The partition function has then the form 



Z n {J) = J dq {ab) dl {ab) dr [ab] dfi [ab] exp [-Nh(q, I, r, fi)] 
where the function h is given by 



(45) 



/^ 2 



h(q,l,r,fx) =- — 



, b | 4 +|r ab | 4 



a (ab) 

Re(V*aa r aa + l* {a b)1(ab) + V*(ab) r (ab)) ~ hi Z(l( ab ) , M[a6] ) 



(46) 



We have extracted the diagonal part in the effective Hamiltonian reminding that q aa = 1 is fixed (this is why we 
didn't include the relative 8- function) . 

The integral (|45ll can be evaluated at the saddle-point. The derivatives with respect to q and r yield 



Substituting these equations into (|46l) . h is 



n/3 2 3/? 2 
h(q,r) = -—+ — 



l(ab) = -2/3 |g( a fc)| q(ab) , 
l-^aa P yaa\ ^aa j 

M(ab) = -2/3 2 |r (afc) | 2 r (ab) . 



]T|r Qa | 4 + 2 5> Qb | 4 + M 

a (ab) 



(47) 



-lnZ(<7 (ao) ,r [a6] ) 



(48) 



To perform the analytic continuation to n — > one has to make an ansatz on the structure of the matrices q and r. 
The free energy is then computed using the equation (|34|l . Using the relation 



lim Z n ~ e -^ mi »['*(9,'-)] 

JV^oo 



and assuming that mxn[h(q, r)] ~ n (i.e. it is small) we have for the free energy 



e -jvmto[h(g,r)] _ j ATmin[/i(g, r)l 

hm ~ hm = mm 

n->0 nA^ n^o nA^ 



lim n h(q, r) — mm[/3(f>(q, r)] 



(49) 



(50) 



Note that the limits n — > and N -^ oo have been exchanged |40|. A reasonable ansatz that we will make is that 
r a b = at the saddle point. Indeed, we are looking for disordered states that are usually characterized by a non-zero 
overlap q and a vanishing magnetization, i.e. the rotational symmetry is not broken. As r is not invariant under 
rotations, we will set r = in the following. Then we have to minimize 



fl2 oo2 

Pct>{<D = - V + T~ E l^l 4 - *"* kZ(g) , 



2n 



a>b 



Z(q)= / e^ a exp 



Re^2/3 2 |g QO | 2 < ? : fc e 



i(v a -v") 



a>b 



(51) 
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for n — > 0. This function is very similar to the one that describes the Ising p-spin glass. 

The vanishing of the derivative with respect to q ao of (3<f)(q) gives the saddle point equation 



q a b 



i{v a -v )\ 



(52) 



where the average is on the measure that defines Z(q). Indeed, performing the derivative with respect to the real and 
imaginary part of q a b, one obtains two equations that can be written in a single saddle point equation for q ao 



3 q ab = (e^-^) + 2. 



q a b 

Itol 2 



R-e q* ab 



i(tp a -lp )v 



(53) 



It is easy to show that Eq. I|52|l is solution of Eq. (|53|) . 

The replica symmetric solution corresponds to q ao = q; in particular q = is a solution of the saddle point 
equations and is the stable one in the high temperature phase. Another solution appear at low temperatures but 
is always unstable (see below). The RS free energy is simply fas — —(3/4 as in the Ising p-spin glass (neglecting 
irrelevant constants). 



B. One step replica symmetry breaking 

The lRSB ansatz is the following: we divide the matrix q a ^ in n/m blocks of side m. The elements in the off-diagonal 
blocks are set to while in the diagonal blocks RS is assumed and q ao = q. The simplest choice is to assume that q is 
real. This is very reasonable due to the rotational symmetry, see e.g. the discussion in |5 lj. p. 894, and moreover in 
this way the constraint q a b = ql a is respected. For instance, we have for n = 6 and m = 3: 



{Qab) 




Then one has 



limn 1 VV h = -(to- l)g 4 , 

n— »0 ^— ' z 



a>b 

and, as the variables ip a in different blocks become uncorrelated 

n/m 



Z(q) 



Udip 

a=l 



a„/3V£ a ^< 



Hf^-f") 



^ a=l 



n/m 



a—1 

where £ is a complex variable and V( = dRc C^ Im C e ~2CC* _ Defining also I — \/2q 3 one has 

n- 1 \TiZ(q) = ~(3 2 q 3 + m- 1 lii f V( f f dipe plKcCe ' V J . 
Introducing the modified Bessel function of the first kind of order 

and noting that it depends on the modulus z = \(\ one finally gets (apart from constant terms) 



nRSB 



r 2 i r°° 

( m , T) = _fL [1 + 3(1 - m)q 4 - Aq 3 ] - - In / Vz J£ 
4 to J 



■(J3lz) , 



(54) 



(55) 



(56) 



(57) 



(58) 



(59) 
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where Vz = ze z l 2 dz. The value of q is determined by the condition d q 4>msB = that gives (see Appendix iBl for 
details): 



J °° vziptfiz) 
q = - 



ijiplz 



Ia(filz) 



^VzI^(filz) 



(60) 



where I\(x) — Iq(x) is the modified Bessel function of order 1. This expression is similar to the Irsb free energy for 
the p-spin model with p = 4, the only difference being the presence of the Bessel functions instead of the hyperbolic 
cosine in the integrals, the domain of integration in z and a z in the integrand. 

The equilibrium value of m is the solution of d m <j)iRSB = 0. At high temperature the solution q — and m = 1 
(paramagnetic state) is the stable one, while for T < T c a new solution with q ^= and m < 1 (spin glass) becomes 
stable. The temperature T c (also called Kauzmann temperature Tk) marks the appearance of the thermodynamic 
glassy phase. 

C. Phase space structure of the model 

Starting from Eq. I|59|l one can repeat the analysis of p3 to derive the full phase space structure of the model at 
the Irsb level. Again, we will reproduce in some details the original derivations for the reader who is not familiar 
with these methods. 

In this class of mean field disordered models, at low temperature, the phase space in disconnected in many metastable 
states, i.e. local minima of the free energy. The number of states of given free energy density / is 0(f) = exp7VS(/). 
The function £(/) vanishes continuously at / = f m i n and drops to zero above / = f max (see e.g. J53|)- The main 
peculiarity of these models is that an exponential number of metastable states is present at low enough temperature. 

One can write the partition function Z , at low enough temperature and for N — *■ oo, in the following way: 



z = e -0NF ( T) „ J- e -PN fa = f fmax df e Jv[E ( /)-/9/] „ e Nmn-pn ; (61) 

Jfmin 

where /* £ [f m im fmax] is such that <&(/) = / — TT,(f) is minimum, i.e. it is the solution of 

-af=T> (62) 

provided that it belongs to the interval [f m i n , f ma x]- Starting from high temperature, one encounters three temperature 
regions: 



For T > Td, the free energy density of the paramagnetic state is smaller than / — TS(/) for any / £ [fmin, f 
so the paramagnetic state dominates (in this region the decomposition (|61|l is meaningless). 



mm 7 J max 1 5 



• For Td> T > T c , a value /* £ [f m in, fmax] is found, such that /* — TS(/*) is equal to f para - This means that 
the paramagnetic state is obtained from the superposition of an exponential number of states of higher individual 
free energy density /*. The phase space is disconnected in this exponential number of regions: however, no 
phase transition happens at Td because of the equality /* — TE(/*) = f para which guarantees that the free 
energy is analytic on crossing Td- 

• For T < T c , the partition function is dominated by the lowest free energy states, /* = fmin-, with Tt(f m in) = 
and F(T) — f m in — TTi(f m i n ) = f m in- At T c a phase transition occurs, corresponding to the 1-step replica 
symmetry breaking transition found in the replica computation. 

In the range of temperatures Td > T > T c , the phase space of the model is disconnected in an exponentially large 
number of states, giving a contribution S(T) = £(/*(T)) to the total entropy of the system. This means that the 
entropy per particle S(T) for Td > T > T c can be written as 

S(T) = E(T) + S mb (T) , (63) 

S V ib(T) being the individual entropy of a state of free energy /*. The task is then to compute the function £(/) at 
fixed T and the equilibrium complexity S(T) = E(/*(T)). 

To this aim, the idea [56j is to consider m copies of the original system, coupled by a small attractive term added 
to the Hamiltonian. The coupling is then switched off after the thermodynamic limit has been taken. For T < Td, 
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the small attractive coupling is enough to constrain the to copies to be in the same state. At low temperatures, the 
partition function of the replicated system is then 

Z m = e -/W*(m,T) _ y^ e -/3JVm/ Q = f ™" ^ gJV[S(/)-/3m/] _ gJV[E(/*)-/3m/*] ^ (g 4 ) 

where now f*{m, T) is such that $(to, /) = mf — TE(/) is minimum and satisfies the equation 

W = T- (65) 

If to is allowed to assume real values by an analytical continuation, the complexity can be computed from the 
knowledge of the function $(771, T) = mf*(m,T) — T£(/*(to,T)). Indeed, it is easy to show that 

r{mX) = ^pn, 

dm 1 (66) 

£(m, T) = £(/*(m, T)) = to 2 * 9 ^ /?$(to,T)] = m ^/ T ) _ ^( T ) _ 

C7TO 

Thus the function £(/) can be reconstructed from the parametric plot of f*(m, T) and £(m, T) by varying m at fixed 
temperature. The equilibrium complexity is simply E(T) = £(m = 1,T). 
Using the replica trick to compute the free energy, 



$( m , T) = -4*og^; = -4 lim ^ = -4 lim , (67) 

one obtains the partition function of nm copies of the system, with the constraint that each block of to replicas has 
to be in the same state, i.e. the replicas must have nonzero overlap. This leads naturally to the Irsb structure for 
the overlap matrix (with to fixed), see Eq. Q54JI . and 

T 
$(to,T) = -— lim[exp [ - f3nmN4> 1RSB (m, q,T)] - l]/n = m4> 1RS B(m,T) . (68) 

iV Tl— HJ 

Note that the hypothesis that the to replicas are in the same state implies that for any value of (to, T) one has to 
substitute in 4>\rsb the nonzero solution of Eq. IjfiOp . q*(m,T). Above T^ this solution disappears as a vanishing 
coupling cannot constrain the replicas to stay close to each other. 
Using Eq.s (|66(l and Ij68(l the complexity as a function of to is 

TS(to, T) = m 2 d m [to-^to, T)} = m 2 d m <p 1RSB {m, q\ T) , (69) 

and the equilibrium complexity is 

3/3 2 (<f) 4 , ^ [°°„_ T im ^, f™-DzI (!3l*z)biIo(pl*z) 



4 J J Q Vz I {pl*z) 



where Z* = y / 2{q*) 3 . 



D. Phase diagram in the [m, T)-plane 



For a given value of T, we can identify four relevant values of to. They are reported in Fig. ^ an d are defined as 
follows: 

1. A solution q*(m, T) 7^ of Eq. 1|6()J) is present for to > m m i n (T). Thus, in the region to > m m i n (T), $(777, T) is 
well defined and we can compute the complexity E(m,T) and the free energy f*(m,T) using Eq.s (|66J) . 

2. The value m*(T) such that E(to,T) = correspond to the solution of the thermodynamics and f*(m k ,T) — 
4>irsb(iti* , T) is the free energy in the spin glass phase. The temperature T c is defined by m*(T c ) = 1. 

3. The function f*{m,T) has a maximum for 771 = md(T). This means that f*(mrf,T) is the maximum possible 
free energy f m ax(T). The states with / = f max are called threshold states |53ll55l |. 
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FIG. 1: (Color online) Phase diagram of the model in the (m,T) plane (see Section IV D in the text). 

4. Finally, one can investigate the stability of the Irsb solution with respect to further steps of replica symmetry 
breaking, following the analysis of |57|. We report the details of the calculation in the Appendix lO It turns out 
that the Irsb solution is stable toward further steps of replica symmetry breaking for m > m s (T). 

The thermodynamic phase diagram of the model is easily deduced from Fig.^ The paramagnetic solution is stable 
for T > T c . The temperature T c is the thermodynamic glass transition temperature, at which a downward jump of 
the specific heat is observed, as in standard first-order transition. Below this temperature the lRSB spin glass solution 
is stable, and remains such down to T = 0, as m*(T) > m s (T) for all T. Thus in this model no Gardner transition 
(transition to a full RSB solution |57|) is observed. 

Some information on the dynamics of the model can also be obtained from Fig. ^ Indeed, at Td (defined by 
rrid{T) = 1) a dynamical transition takes place [53|: correlation functions are expected to develop an infinitely long 
plateau and the system becomes dynamically trapped in one of the exponentially large number of states that appears 
at Td, as discussed above. 

Between the static transition temperature T c and the dynamic one Td the phase space has a non-trivial shape: it 
is disconnected in an exponential number of states Af(T) = exp7VE(T), where E is the configuration entropy (or 
complexity) of the system. In Fig. |2 the quantity E is reported as a function of T. The point T c at which the 
complexity vanishes E(T C ) = signals the appearance of the thermodynamic phase transition. 

However, it is not clear what is the dynamic behavior of the model if quenched from T = oo to T < Td- Indeed, 
as already observed in the Ising p-spin glass model, the threshold states always lie into the lRSB-unstable region, i.e. 
md < m s for T < Td- This means that the Irsb ansatz is unable to give the correct prediction for these states below 
Td and one need to consider further steps of replica symmetry breaking. The investigation of the dynamics of the 
model after a quench below Td will be the object of future investigation. 

V. ROUTE TO THE EXPERIMENTAL OBSERVATION OF THE GLASSY TRANSITION 



The investigated model exhibits a dynamical transition at Td and a thermodynamic phase transition at a lower 
temperature T c (characterized by one-step replica symmetry breaking scenario), very similar to p-spin glass models. 
Let's turn our attention to the physical interpretation of these transitions. We recall that the "temperature" T 
introduced in the model is defined as the ratio between the bath temperature Tb at h, measuring the optical noise in the 
system, and the pumping rate V, so we can vary the latter to explore the phase diagram of the system. 

Starting from low pump intensity (high temperature) and increasing V (decreasing T), different interesting phe- 
nomena take place. The relevant quantities to look at in experiments are correlation functions in time domain, i.e. 
the self correlation functions of a specific frequency (u> m ) component of the electric field in the cavity (for example 
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FIG. 2: The equilibrium complexity £(T) as a function of the temperature. 



via heterodyne experiments, see below): 

C(t,u> m ) = (a m (t + T)a* m (T)) T = A 2 m (exp{i[<^ m (i + r) - ¥>m(i")]}) T 



(71) 



where the (. ..) T is the average over the time origin r. At high T, because of the fast dynamics of the phases, 
C(t,u> m ) decays to zero on short times. On lowering T the dynamics of phase variable f m (t) becomes slower and 
slower and C(t,uj m ) is expected to decay towards zero in longer and longer times. At the dynamic transition point 
Td, the dynamics of the ip's becomes non-ergodic, they are no longer able to explore the whole phase space and 
the function C(t,u) m ) decays towards a plateau: the mode's phases tp m (t) are locked to some fixed random values 
("random mode-locking") and oscillate around these values. The fact that the complexity £ is different from zero 
at the dynamic transition point implies that there is an exponentially large number of possible values for the locked 
phases and, correspondingly, many different time structures of the electric field in the random laser. 

More technically, the region between Td and T c is dynamically not-accessible, due to the mean-field character of 
the interactions. Only for short range models, where activated processes become possible, this interesting region can 
be explored. We can expect that, if the mean-field approximation is not fully verified by the random lasers, see the 
discussion in section IIII Al the system would be able to enter in the activated (or Vogel-Fulcher) regime (to use a 
terminology familiar in the liquid-glassy physics) and the correlation functions decay to zero at long times after the 
plateau. In this region the relaxation time is found to scale as r = To exp[D / (T — T c )] and only at the thermodynamic 
phase transition point T c the system is really locked in the ideal glassy state (the relaxation time becomes infinite). 
Similarly to what happens in p-spins and structural glasses, interesting phenomena as aging, memory effects, and 
history dependent responses are expected to take place for T c < T < Td, see e.g. |49|, |58|, |5S| for recent reviews. 

It is worth to remark that in the region where these phenomena are expected to happen, the relaxation time of the 
system is larger, by many orders of magnitude, than the typical microscopic time scales. For instance, in molecular 
glasses where the typical time scale is To ~ 10~ 12 s, the relaxation time of the system can be as large as 100s already for 
T g ~ 1.5T C . These systems cannot be equilibrated close to T c because of the exponential divergence of the relaxation 
time for T — ► T c . In random lasers, the typical time scale for photon dynamics is To ~ 10 _14 s, so one will gain orders 
of magnitude in time. This should allow to equilibrate the system closer to T c . Even if the decrease in T — T c that 
one can achieve in this way will not be substantial, it could be en oug h to test the predictions of some recent theories 
of the glass transition in presence of short range interactions |43, |5g ■ 

Moreover, we recall that, as discussed in section IlII Al in random lasers it is possible, in principle, to tune the 
interaction range (by tuning the localization length), and thus the relevance of the activated processes. This is similar 
to what has been done theoretically by considering the Kac limit for spin glasses [6(| and might allow to explore the 
crossover between the mean field and the activated regime and to shed light on some debated issues concerning the 
nature of the spin glass phase in real systems. Before concluding, it is important to provide an order of magnitude 
estimate of physical quantities involved in the experiments and to address some possible experimental frameworks. 
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A. Order of magnitudes 

For the sake of concreteness we will focus, as an example, on recent experiments in random lasers realized by 
scatterers (e.g. zinc oxide powder dispersed in a solvent doped with a dye, as e.g. in |2l|). These experiments 
employed pumped source beams (see below), for the moment we show that the glassy transition can be expected for 
the currently employed pump power levels. We stress once again that our model is sufficiently general to embrace a 
much wider variety of disordered amplifying systems. 

We start from the coupling coefficients, which are given by Eq. (j25(l : the fields are real valued numbers, their 
modulus scale as \E\ = 2V^ 1 / 2 / \Je n^ due to the normalization (n is an average refractive index); their sign can be 
"embedded" in the sign of the x coefficients; additionally uj s oc cj = 2irc/\. Hence, omitting indexes, (g) — and 

8 "° 2 ' x(r)dV (72) 



v e o n o Jv 
We need {g 2 }, and we assume 

( X (r)x(r')>=X^(r-r') (73) 

with xo a typical nonlinear susceptibility value, L r a characteristic length for the disorder and S the coarse grained 
Dirac delta. Thus 

<-' 2 > * (s4?) 2{ L- x{r,dv 'L x{ '' )dr) - (Mz) xlLlv (74) 



We have then for the standard deviation 



Suj 2 xoLl /2 I 



For the sake of simplicity, let us consider a box with volume V: the number of modes per unit of volume and unit 
of frequency is given by (this is just an estimate, as in nanostructured optical cavities the density of modes can be 
enhanced or depressed with respect to a standard box Q) 

PW) = ^ (76) 

hence for TV modes in a wavelength range AA (A = cjv) it is 



Eq. (|77|l is used in l(75)) and gives 



N^-^V (77) 



^ , S^ 3 / 2 ^ 2 ^ 2 (AA) 3 / 2 1 

V{g } " 44 a^ a^' (78) 



which scales as N 3 ' 2 as anticipated. Next we have to consider the coefficients G sp qr — gspqrA s A p A q A r 

i^V/^xoLl^n 1 ^ 2 (AA) 3 / 2 (A 2 ) 2 



V(G}- ^ A^A^' { ' 

and the coefficients J = G/(go(A 2 } 2 ) with variances 8/N 3 . We can hence determine go as 

„ Sn^nl^xoLr (AA) 3 / 2 

90 = 1 2 ^~ (80) 

and, finally, the adimensional (3 — \/T is given by 

^ulxoLl' 2 (AA) 3 / 2 (A 2 ) 2 

£n A° k B 1 bath 
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Remembering that the average energy per mode is cuo(A 2 ), we obtain, at the transition, the threshold value (denoting 
Trf ~ 0.435 the normalized temperature at the dynamic transition, see Fig. |2J): 



£ RSB = MA 2 ) = / *BWgA° 

The noise temperature can be taken as due to the spontaneous emission, which is typically the dominant contribution: 
following 42], 2k B T bath = h(N 2 /{N 2 - iVi)) t /r = fi/r, with r the average life time per mode, and (N 2 /(N 2 - N x )) t 
the population inversion at lasing threshold. Taking typical values from the reported experiments, (AA = 100 nm, 
A = 630 nm, n = 2,L r = 10 nm, r = 100 fs) and for the susceptibility xo = 10~ 27 CmV" 3 H3 it is £ RSB = 1CT 16 J. 
Assuming a pumping beam with peak power Prsb — NSrsb/t, gives Prsb — 0.1 W with N = 100, which focused on 
the typical area of 100 fim 2 provides the typical values for the peak pump intensities used in the experiments (~ 100 
kW/cm 2 ). Thus we expect that the "glass transition" can be observed within the currently available experimental 
framework. 

B. Continuous- wave random lasers 

All the theory reported in this manuscript makes reference to continuous- wave (CW) RLs, hence we start discussing 
this kind of systems. Experimental investigations of CW-RL were already reported in nano-powders (see |6l| and 
references therein); alternative experimental geometries include disordered photonic crystals |2|, as membranes or 
multi-layered systems with gain provided e.g. by quantum wells in semiconductor materials. In these integrated high- 
index contrast geometries multi-mode random-cavity lasers can in principle operate in CW regime, and this opens 
the way to a comprehensive experimental analysis of the dynamics of the laser emission. The experimental setup 
can follow previous investigations of the noise figure of standard semiconductor lasers (see e.g. 62]). RL emission 
is collected, and filtered in a narrow band in order to select one or few modes. At the mode locking transition the 
intensity signal I(t) is expected to switch from a random noise superimposed to a CW value, to a largely modulated 
line-shape displaying a random sequence of disordered pulses. Heterodyne measurements are employed to extract 
phase and amplitude noise from which the dynamics of the amplitudes of the modes and their phase are extracted, 
as well as the coherence function C(t,ui). As discussed above, the fact that amplitude fluctuations make a negligible 
contribution to the field autocorrelation is well known from laser theory [62|, hence C(t, u) gives information on the 
phase-dynamics. 

Specifically, before the glassy transition (low pumping rate) the line-width of the laser modes is wide and the 
autocorrelation C(t, us) of the mode signal is expected to decay with a single exponential time constant, corresponding 
to a Lorentzian line-shape of the noise spectrum of the field. At the glass transition the laser dynamics slows down, 
and this result into a slower decaying of C(t,oj), with the appearance of multiple time scales, and eventually to an 
ergodicity breaking corresponding to a plateau in the C(t,u>) signal (as sketched figure^)- 

Further information is retrieved by phase demodulation (e.g. by employing the usual combination of a limiter and 
a discriminator |62|) whose output is the instantaneous frequency d<j) m /dt, whose power density spectrum can be 
determined by a spectrum analyzer. When the phase of the filtered modes are locked, they oscillate around one of 
the many equilibrium values. Correspondingly, the spectrum of the phase noise for each mode is expected to switch 
from a wide line to a narrow one displaying modulation side bands, which are due to the fact that the phases display 
small oscillations around one of the many phase-locked states (as sketched figure [5J}). 

C. Pulsed random lasers and speckle patterns 

Most of the reported experiments on RLs have been done by using pumped laser beams, with pulse duration from 
tens of picoseconds to tens of nanoseconds. One could argue if the mentioned mode-locking transition can be actually 
observed in these regimes. 

First of all we observe that our theory deals with the phase-dynamics of the (quasi-)modes of a RL. As far as the 
laser reaches a steady state for the amplitudes (and this is expected to happen in the leading edge of a nano-second 
pump pulse, taking into account typical lifetimes; see e.g. |20| for an extensive discussion of the mode amplitude 
dynamics), the phases are expected to vary on the To ~ 10 fs time-scale. Indeed, in the framework of the Lamb's 
two level theory, they are affected by the dynamics of the resonant medium polarization, i.e. by the time-constant of 
the off-diagonal density matrix elements, whose inverse is the "dipole dephasing rate" which is around 10 14 s _1 (see 
e.g. |34()- The latter time scale is the "elementary" time scale for the dynamics of the phases, that corresponds in 
molecular systems to the typical time scale of atomic vibrations tq ~ 10~ 12 s. 
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FIG. 3: Sketch of two experimental signatures of the onset of a glassy transition of light in random lasers: (a) The mode 
coherence function develops a plateu at high pumping rates, denoting ergodicity breaking, (b) The power density spectrum of 
the mode instantaneous frequency displays modulation sidebands; this corresponds to the transition from random-mode-phases 
to phase-locking in one the meta-stable states; the side-bands are due to small oscillations in these minima that result into 
frequency shifts of the mode resonances. 



In a first approximation, the arrival of a pulse produces a fast variation of the "temperature" from T ~ oo (before 
the pulse) to a value of T < Td (subsequently after the arrival of the pulse). This corresponds, in the language of 
molecular glasses, to an instantaneous (i.e. on the scale of To) quench of the system from infinite (or very high) 
temperature to below Td- On a very general ground, it is expected that the system will then start to age |49ll59j. in 
the sense that the relaxation time r of C(t) will increase with the time t w elapsed after the arrival of the pulse (the 
quench). For systems that are in the class of the p-spin model it is generally found that, for t w ^> To, r{t w ) ~ t w (see 
as a striking example Fig. 5 in [63]). This means that for a pulse duration t w ^$> tq ~ 10 fs, the relaxation time of the 
phases will be of the same order of the pulse duration and they will appear to be frozen on this time scale. Hence the 
mode-locking process should be observable in standard random laser for ns pump pulses. 

Note that, incidentally, the validity of the CW approximation has been thoroughly discussed in a recent paper for 
nano-second pulsed laser |24j . 

Additionally, the pulsed regime favors the investigation of the correlation between different laser shots with fixed 
disorder. In the presence of many modes, due to mode interferences and complex phase modulations, the peaks in the 
spectrum are expected to largely vary from pulse to pulse. This is also due to the fact that with a pulsed pump, in 
a thermodynamic language, the system is first "cooled" (i.e. the average energy per mode is increased as the pump 
power increase) and then "heated" (as in the trailing edge of the pump pulse) J7Q| . Between two subsequent pulses 
the system is at infinite temperature and will rapidly loose memory of the previous state: as a result, in the presence 
of an exponentially large number of thermodynamically equivalent states, the system settles in a different minima 
from pulse to pulse and correspondingly the phases will largely vary from shot to shot. During the laser oscillation, 
the phase-modulation corresponding to the mode-locking process should be visible. Similar phenomena (namely a 
large variation of the emitted spectrum from pulse to pulse) were already reported in the literature (as in Ref. |21|'). 

A note on the speckle pattern of the emitted light is in order. Indeed the speckle pattern is determined by the 
phase difference between the modes, hence the spatial distribution of the emitted light is expected to largely vary from 
one shot to another of the random laser, when the pump power is above the threshold for the glass transition. It is 
important to emphasize that other authors predicted an exponentially large number of speckle patterns in nonlinear 
random media |64|. however, in that case, the leading mechanism was the non-resonant Kerr effect (incidentally, our 
model Eas. (|31|) also applies for Kerr media as will be discussed elsewhere). 



VI. CONCLUSIONS 



We derived a statistical model for the mode phases in a random laser. The obtained Hamiltonian resembles that 
of some spin glass model (the p-spin model with p = 4 |55l 15 7L 65, 66] ), and can also be thought as a generalization 
to the disordered case of a toy model Hamiltonian recently studied (the ^-trigonometric model |67|1. The relevant 
parameter for exploring the phase diagram is a scaled temperature T, the ratio between the "true" bath temperature 
and the square pumping rate (the energy stored on average in each mode). Using standard statistical mechanics 



18 

techniques of disordered systems (i.e. replica method), we predict the existence of a dynamic transition at Td and 
of a thermodynamic phase transition at T c , characterized by one-step symmetry breaking scenario. Between the 
two temperatures, the appearance of an exponentially large number of states is expected. This corresponds to the 
existence of a random-mode locking transition in random lasers: looking at self-correlation functions of a specific 
frequency component of the electric field in the cavity, one should observe a non-ergodic behavior at Td, i.e. the decay 
towards a plateau, where the phases are locked at random values. The mode-locking will happen in a configuration 
of the phases that depends on the history for a given sample, because the system will reach a different metastable 
state depending on the initial state, with large sample-to-sample fluctuations. Thi can be observed by looking at the 
frequency fluctuation spectrum or at the speckle pattern of the emitted light. 

The logarithm of the number of these possible random configurations of the phases is given by the complexity (see 
Fig. [5Jl times the number of active modes. As long as the physical realization of the random laser is well described by 
the mean-field Hamiltonian, the system is not dynamically able to explore the phase space for temperature T < Td 
(or pumping rate higher than that corresponding to the dynamic transition) and remains trapped always in the state 
reached at Td- However, taking into account the fact that the mean-field scenario could be a too "crude" approximation 
for real random lasers, we expect that the system will be able to explore on a long time scale T c < T < Td region, 
where aging, memory effects, and history dependent responses are expected to take place. 

The expert reader in nonlinear optics or Bose-Einstein condensation will certainly recognize in our model a typical 
system for many-modes interaction processes (e.g. solitons, parametric processes, supercontinuum generation ...); we 
believe indeed that our result are also relevant in many branches of modern nonlinear physics involving disordered 
systems. Spin glasses have been defined has the "most complex kind of condensed state" J68J, we are convinced that 
there is no difficulty in accepting the emission of random lasers as the "most complex kind of light" . 
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APPENDIX A: EQUATIONS FOR THE AMPLITUDES AND FOR THE PHASES 

Here we will show that the Hamiltonian (|3U[) can be derived also directly from the equation of motion for the 
complex amplitudes a s , Eq. I|26(l . by assuming that the dynamics of the phases is much faster than that of the A s . 
After l|2ti[) . using g sp qr = g^pqr + wlpqr an< ^ Vs = V^ + *%> the equations for the amplitudes are: 

-jf = ~ 2 Yl A p A Q A r [dspqr cos(ip 8 + tp p - ip q - ip r ) + g\ pqr sin((/p s + ip p -<p q - <p r )] + 

pqr V / 

+ (7 S -a 8 )A s + rjf- cos(ip s ) + t][ sm{ip s ) 

By assuming the phases as rapidly varying with respect to the amplitudes, these can be averaged out and the "free 
run approximation" equations |20l. l26t l28j are retrieved, providing the average energy in each mode £ s : 

AC „R „R I n R 

° s o/ \ c "ssss c2 i c \ .Jsrsr ' i* ssrr c / \n\ 

—TT = 2(7 S - a s )t s t s + t s > t r ■ (A2) 

at uj s ^— ' uj r 



Note that Eqs. (|A1|) take into account random cross- and self-saturation effects, by the terms weighted by gf rsr +gf srr 
and gf sss respectively, as well as the fact that the decay rates and gains are expected to be different for each mode. 
We model this circumstances in the text by taking the amplitude dependent G-coefficients as gaussianly distributed. 
Next, we find the ruling equation for the phases after l|2()|l : 

A s — t^- = --^ApA q A r {gl pqr cos{(p s + <pp-(p q -<p r )-gi), qr sm((p s + <p p -(p q -(p r )}+rilco8(<p s )-TiK siu(ip s ) . (A3) 



pqr 



Denoting £ = uj (A 2 ) the average energy per mode, we multiply 1A4() by ujqA s and take ojqAJ. = £ (while being 
lo s = ujq, as outlined above). 
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The equation for the phases l|A3() are equivalent to the following 



dip s 
dt 



gff, 

dip 



V i „(¥>) 
'/s 



with 77s = [ryf cos(< J 5 s ) - ryf sin(y> s )]/A s and 



*.-£!: 



1 >r-^ uj Q A s A p A q A r j 



8 



[g shx{cp 8 + <p p - (f q - ip r ) + g cos((^ s + ip p - ip q - <p r j\ 



(A4) 



(A5) 



spqr 



where we exploited the symmetries for the real part which are the same as discussed above, and where the imaginary 
part satisfy g J apqr = -g{ rap . 

Eq. (|A5Jl is cast in the form H v = luqH/2£o with 



H = - Y^ A s A p A q A r [gR pqr cos{<p s + <p p - ip q - ip r ) - gl pqr sin(^ s +(p p -<p q - <p r j\ = -Re 

spqr 



Q spqr&qQ'rQ'pQ' s 



spqr 



(A6) 



Eq. (|A4|I is a Langevin system for the phases, and being (due to the fact that the noise r\ is assumed to vary 

on a much faster scale than the phases ip) (r] p (t)r]q' (£')) = i^okBTbathS pq S(t — t')/£o, its invariant measure is 
exp(— 2H v £ /ujQkBT} )at i l ) — exp(— H j feeT;,^) j which is identical to the previous one. In the case of a complex 
coupling the Hamiltonian will include also the second term in (|A6|) . The replica analysis of the generalized model is a 
generalization of the one concerning Eq. 1|27(1 (see section IIV() and provides the same outcome for what concerns the 
existence of a RSB transition. 



APPENDIX B: SELF-CONSISTENCY EQUATION OF 1RSB SOLUTION 

Here we derive the self-consistency equation for q, Eq. (|60() . It is obtained imposing that the derivative of the free 
energy in Eq. I|59[) with respect to q vanishes, d q 0<fii as.B(m, q) — 0: 



30 2 q 2 [(l-m)q-l] 



1 J^Vzd^iPlz) 



Now, d q I™ = mI™~ l d q Io and d q Io = (d a Io)(d q a) — Iid q a, where a = (3lz — (3^/2q 3 / 2 z. Then we can write: 

\/2 Jo 



Vz W = ^ 



Vz zl^- l h = ^=- I dze- : 'l 2 dzizl™- 1 !,) 



(Bl) 



(B2) 



having used the identity ze z ' 2 — —d z e z ' 2 and integrating by part. Performing the derivatives and using d z = 
(d z a)d a and the property of Bessel functions d a I\ = Iq — I\/a, we have: 



m- 1 / Vz d q I r n = 3/3 2 g 2 / Vz 7 m 1 + (m - 1)4. 
Jo Jo l A) . 

Substituting in Eq. I|B1(I we then obtain the self-consistency equation: 

J °°VzW(/3lz 



(B3) 



I {plz) 



J^VzI^iPlz) 



(B4) 



APPENDIX C: STABILITY OF 1RSB SOLUTION 



In this Appendix we discuss the stability of the Irsb solution. First we have to compute the Hessian of (f>(q) 
evaluated in a solution that verifies the saddle point equations q a b = (cos(<y3 a — ip )), where we assume that q a b is real. 
This assumption is motivated by the analysis of [54]], p. 894, where it is shown for the case of a two spin interaction 
that even the full RSB solution verifies q a b real for all ab. Still this is an assumption as in principle there could be 
solutions such that q a b has an imaginary part for a ^ b, see again |54j . 
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Considering a perturbation 6q a b around the solution, one has, differentiating Eq. (|51|l : 
d 2 /3<j> 6/3 2 



G a b,cd 



dq a bdq c d 



-q 2 ab [5ab,cd - &[3 2 q 2 cd [{cos(<p a - ip b ) cos(ip c - </)) - qabqcd] 



(CI) 



The condition q ab = q£ a implies 5q ab = Sqb a ■ The matrix G is symmetric under the exchanges a <-» b, c <-> d. When the 
matrices G is evaluated in the Irsb solution, it is easy to see that it becomes a block matrix that has non-vanishing 
elements only if (ab) and (cd) belong to the same diagonal block related to one of the diagonal blocks of the matrix 
q a b- Thus we can restrict to consider perturbations of one single block. With this restriction, substituting the lRSB 
structure of q a b, and neglecting the irrelevant prefactor 6f3 2 n~ 1 q 2 , G has the following elements: 



P = G abMb = 1 - 3/3V 

Q = G a b,ad = —3/3 q 

R s G ab , cd = -6/3V 



jvchmc\) m ' 2 hmc\f _ „ 2 

JV(I ([3l\C\) m \ 

jvcio(pi\c\) m - 3 i2mc\)hmc\) 2 

q Jvcio(pi\C\) m 

fDCIo(Pl\C\) m - i h(Pl\C\) i 



2q 2 



(C2) 



fv(i (i3i\c\y 



Following the analysis of [fijj, the relevant eigenvalue of G that eventually becomes unstable is L = P — 2Q + R, i.e. 
the stability condition is 



/*>CWICI) 



1 - 6/3 V 



m ) 1 
2 



1 



MffllCI) 2 



MffllCI) Ji(ffllCI) 2 

Ia(f3l\C\) Io(MC\r 2 



Jv(i (f3i\C\Y 



>0 



(C3) 



that is 



/2?Cio(/3i|CI) 



6f3 2 q 



2„2 



> 



mjl 



MfflKi: 



M/3/ICI) 2 



/ (/3;|ci) /o(/3/|CI) 2 



fD{Io(fil\C\) r 



(C4) 



Once the solution of the saddle point Eq. (|B4|I is substituded in the expression above, one obtains the condition 
m > m s (T), see Fig. ^ 
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